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Given a natural number N , one may ask what configuration of N points on the two-sphere mini¬ 
mizes the discrete generalized Coulomb energy. If one applies a gradient-based numerical optimiza¬ 
tion to this problem, one encounters many configurations that are stable but not globally minimal. 
This led the authors of this manuscript to the question, how many stable configurations are there? In 
this manuscript we report methods for identifying and counting observed stable configurations, and 
estimating the actual number of stable configurations. These estimates indicate that for N approach¬ 
ing two hundred, there are at least tens of thousands of stable configurations. 
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I. INTRODUCTION 

“ Computer trials indicate that in the range 70 < A < 112, the number of distinct configurations associ¬ 
ated with each value of N grows exponentially, i.e., M(N) = 0.382 exp(0.0497A). If this trend is sustained 
for larger values ofN, identifying global minima among a large set of nearly degenerate states for complex 
systems of this type will pose formidable technical challenges .” T. Erber and G. Hockney™ 95 

For a natural number A, we denote by qjn - {ri,..., r^/} any configuration of A distinct points on S 2 . 
For a non-negative number s, one can ask what configuration minimizes the energy 

( r~ s when s > 0 

E s ((l>n) := 2_J Zj ks ^ Ti - r ^’ where ks ^ = | 

1=1 j=i+ 1 [ -logr when s = 0. 

For s = 1 this is known as the Thomson Problem Tho04 . At first glance this problem seems remarkably 
simple, yet there is not a simple solution. In fact, Smale has identified a variant of this problem as worthy 
of focus for this centur^ 2 ^. 

The current theoretical progress is limited. It is known that for any A and s a globally minimal con¬ 
figuration exists. For a few special cases of A and s there are rigorous proofs that certain configurations 
are globally minimal. Finally, there are some asymptotic estimates for the minimal energy as a function of 
A. In this last category Polya and Szego p ®^, using measure-theoretic arguments, established some elegant 
estimates when s is less than the dimension of the set on which the problem is posed, e.g., 2 in the case 
of S 2 (also cf. Fandkoff^S (pp. 160-162)). Hardin and Saf t HS05 !, and Borodachov, Hardin and Saf$ BHS08 
established similar results when s is greater than or equal to the dimension of the set in question. 

While theoretical progress is difficult, the analyticity of k s makes this an inviting problem for numerical 
optimization, particularly gradient-based optimizations. Work in this area goes back to at least 1977 MK S23 
and there are two efforts that particularly motivated this current effort. The first is Erber’s and Hockney’s 
reports and commentary™ 91 ™ 95 ™ 97 on computational experiments for the Thomson problem for A up 
to 65 and A up to 112, where they provide some initial estimate for the growth in the number of stable 
configurations as a function of A. This work also includes estimates for the first two terms in the asymp¬ 
totic expansion for the minimal energy. The first term is in agreement with the earlier work of Polya and 
Szego, and the second term was later identified in a formal conjecture by Kuijlaars and Saff KS98 and later 
generalized to a large class of two-manifolds CGS+13 . The second effort is the work by Wales and Ulker wu ®, 
and the work by Wales, McKay and Altschulei WMA09 that led to the Cambridge Cluster Database, which 
reports, for many A and s = 1, the lowest known energy for the Thomson problem. 

A significant challenge for numerical optimization is that many configurations are locally minimal, i.e. 
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stable, with respect to E s , but not globally minimal. This motivated us to attempt to answer the question: 
how many stable configurations for a given N and s are there? An earlier work that answers a similar 
question is that of Hoare and Mclnnis, who identify the distinct stable clusters of modest numbers of point- 
particles interacting through Lennard-Jones and Morse potentials™ 76 . For the Lennard-Jones potential the 
number of stable clusters grew rapidly. The present work estimates this growth for the generalized Thomson 
problem for N up to 180 and for s = 0, 1, 2, and 3, and reports some methods we found useful. 

A central question is whether the number of distinct local minima within the energy landscape grows 
exponentially with the number of points. Stillinger and Weber present the following informal argument 
suggesting an affirmative answei SW82 (p. 980). If one can convert one stable configuration into another 
with changes that are localized in space to within a fixed number nearest neighbor lengths, then, as the 
number of points grows, so should the number of available independent changes, and the number of stable 
configurations grows will grow exponentially. If the number of stable configurations does not increase 
exponentially with N , then this would suggest that changes from one stable configuration to another cannot 
be accomplished with only localized changes. 

Our work began by generating a large library of stable configurations, and we describe our methods in 
Section [IT| In doing this, we found that our optimization program would find rotations and reflections of 
the same stable configuration. In response, we used graph-isomorphisms of the Delaunay Triangulation as 
a means to recognize quickly a particular stable configuration. This method accelerated our work consid¬ 
erably, but there are some subtle ways it can fail. In particular we found two distinct stable configurations 
whose Delaunay Triangulations share the same graph structure. This is described in Section [Dl| Even after 
months of numerical experiments running on many compute cores, the fraction of the most recent experi¬ 
ments that generated new configurations never dropped to zero, making it clear that there are many stable 
configurations we did not see. The problem of estimating the number of configurations we didn’t see is 
an example of the broader “unseen species” problem, which arises in many settings such as linguistics and 
ecology. We apply a method developed partly by linguists to provide estimates in Section [TV] for the actual 
number of stable configurations. 


II. STABLE CONFIGURATIONS 

A. Optimization 


We used an iterative unconstrained optimization strategy described in Section III.A of our prior 
work CGS+13 to generate candidate stable configurations. The method consists of non-linear conjugate gra- 
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dient (NLCG) with line minimization and, when that method no longer made progress, Newton’s Method. 
Our experience is that NLCG with line minimization was most effective up until near the end of the op¬ 
timization. Near the end of the optimization calculation, presumably when the configuration is at a point 
where the objective function E s is locally quadratic, Newton’s Method would often make progress when 
NLCG could not. 

When computing the energy, E s (<jl>n), one has roughly N 2 /2 summands that vary widely in range, and 
a direct summation can lead to roundoff errors. We controlled for this error by logarithmically binning our 
summands and only adding the content from the same bin. This allowed us to ensure that we never added 
two numbers whose ratio was more than two or less than one half, until the end when we summed the 
contents of the bins from lowest to highest. Because we could bound the error for summation in a given 
bin, and because we could count the number of summations in the bin, we were able to estimate the error 
in our sums. This approach follows the work of Highairl—and Demel and Hida DH03 . 

B. Testing for Stability 

Given a candidate stable configuration, we use the criteria described in Section III.B of a previous 
publication CGS+13 to test for stability. The central assumption in this criteria is that our iterative optimization 
strategy will produce a candidate configuration a> c N that is close enough to an actual stable configuration 
so that the gradient at which is zero, may be expressed as a linear expansion of the gradient about qj c n . 
That is 


0 - VE s (ojn) « VE s (cj%) + V 2 E s (oj c n )(o) N - oj c n ), (1) 

where VE S and V 2 E S are the gradient and the Hessian respectively of the objective function with respect 
to the 2N angular free parameters. If this approximation were exact, it would allow us to bound the term 
a)N ~ w c N , where subtraction is applied to the 2ALdimensional space of configurations. Conceptually the 
calculation is 


V E s (OJ]y) VE s ((0]y) — ((On £U/y), 

11V E s (0J N ) 1^ ||V£' i y(6L>^y)||2 > ||^iV — — IK — KL * 

Here || • Ih is the unnormalized two-norm allowing the bound of the infinity-norm. 

The Hessian is not invertible, however. For our choice of coordinates there are three rotations of the 
sphere that do not change the relative distance between the points and hence don’t change the energy. 
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While there are choices of coordinates free of such rigid rotations, those coordinate systems degraded the 
performance of NLCG. Consequently the three lowest eigenvalues of the Hessian are zero. The gradient 
has no projections along the corresponding eigenvectors, and we may choose a rotation of so that the 
difference con ~ ^> C N similarly does not project along these eigenvectors. We let TA n denote the fourth lowest 
eigenvector of the Hessian and then we have 


l | V ^(^)|| 2 


> - M %\\ 2 ^ ||^N - 


Nil 


Change in angle on the sphere bounds from above change in position, and so 6)# - ^nWoo P rov ides a bound 
on the distance between corresponding points in the configurations and cj c n . 

Our criteria for stability is that 


|| V £,( o^)|| 2 ^ min r ,^ r . e ^ |r, - r ; | 

“ 10.000 ’ 

which, in conjunction with the assumption that error in the approximation in Eq. ([!]) is negligible, leads to 
the conclusion that no point in a> c N is further from the corresponding point in a) N by more than one ten- 
thousandth the minimum pairwise separation of the points in cj c n . An important consequence of this is that 
if two configurations satisfy Eq. ([2]), and if there is a rotation and reflection that aligns them to within one 
five-thousandth of both of their minimum pairwise distances, then, we say that they are instances of the 
same stable configuration. As previously noted CGS+13 this criteria relies on bounding the infinity-norm with 
the unnormalized 2-norm. Such a bound is tight only when all the components except one are zero. The 
implication is that the maximum difference between a point in our candidate configuration and a true stable 
configuration is likely considerably less than one ten-thousandth of the minimum pairwise separation of the 
points within the configuration in question. For candidate configurations that we believed were instances of 
the same stable configuration, we could often align them to greater accuracy. 

The condition in Eq ([2]) is a useful, reasonably motivated, heuristic for marking a configuration as stable. 
More rigorous bounds would require estimating the error in Eq. ([!]). 


III. DELAUNAY TRIANGULATIONS AND GRAPH ISOMORPHISMS 

Since the energy E s depends only on the distances between points, it is invariant under isometry. How¬ 
ever, that two configurations of points have the same energy does not ensure that there is an isometry 
between the two configurations. With this in mind we only called two configurations the same if we could 
find an isometry that mapped one configuration onto the other to within the tolerances described in the 
previous section. 
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This leaves the question of how to search for an isometry between two configurations of similar energy, 
which we’ll denote here as co l N - {si,..., s^} and cl> 2 n - {ri,..., Y 2 }. A simple approach is to apply Al¬ 
gorithm [I] described in this manuscript. While there are some optimizations such as, at line 5, first testing 
that |r/ - ry| = |si - S 2 I, this algorithm is expensive and must be applied to every pair of configurations with 
similar energy. 
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ALGORITHM 1. A simple way to search for isometries between coL and oj 2 


Isometry Found False. 


2 min 


1 \ Si -Sj\ min 2 |r~r ; | 


10,000 


10,000 


for r, € (o z N do 

for Vj e oj 2 N \{ri] do 

if there is a rotation of (jj 2 n so that r ; = Si e co l N and so that r y - = S 2 e (o l N to within s then 
if this rotation is such that ||<uL - cj 2 \\ < s then 

m jy w 1 loo 

Isometry Found <— True, 
else 




the reflection of the rotation of u 2 about the plane defined by si, S 2 and 0. 


if < - a^ll 


< s then 
Isometry Found <— Tme. 

end if 
end if 
end if 
end for 
end for 


A. Delaunay Triangulations 

We found a more effective algorithm was to look for isomorphisms between the graphs formed from the 
extremal edges in the Delaunay Triangulations of the configurations in question. For brevity we shall refer 
to this as an extremal triangulation. Essentially we are looking at the edges in the set of triangles that make 
up the surface of the smallest polyhedron containing a configuration, ojn. 

To compute the extremal triangulation we used the QHULL software package BDH96 . One immediate 
observation was that certain configurations did not have unique extremal triangulations, for example, the 
configuration with the lowest observed energy for N = 24 and s = l, shown in Figure [I] The four points 
displayed toward the middle of the image are the vertices of a square, and either diagonal can be part of a 
valid extremal triangulation. Because there can be degenerate extremal triangulations, the assumption that 
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FIG. 1. This is one of many possible extremal triangulations for this configuration of 24 points. 

distinct extremal triangulations indicate non-isometric configurations is not, in general, correct. A simple 
test for non-degeneracy is to compute the set of unit normal vectors for the extremal faces, and to make sure 
that the dot-product of any two is bounded away from one. 

In the case that a configuration has a non-degenerate extremal triangulation, the edges of the triangulation 
and the points in the configuration form a graph, and this graph is invariant under rotation and reflection 
of the underlying configuration. In the degenerate case, a rotation or reflection may lead QHULL, due to 
round-off errors, to find a different, but equally valid, extremal triangulation. 

B. Graph Isomorphisms 

A graph on a sphere is a planar graph in that, by choosing one face to be mapped to the unbounded 
component of the plane, the graph can be mapped onto the plane. In doing this the edges that bound this 
face are retained, and the graph structure is preserved whether the graph is embedded on § 2 or R 2 . There 
are efficient algorithms to determine isomorphisms of planar graphs and we use one following the work of 
Lin j Lin 8 Q l. The approach is to generate a tag for each graph with the property that two graphs are isomorphic 
if, and only if, the two tags are the same. The cost for finding isomorphisms between M instances of graphs, 
or for finding isometries between M configurations, can be written as 

CiM + C 2 M 2 . 

In our approach C 2 is the cost of a searching for matching tags, i.e. string comparisons, while C\ is the cost 
of generating the tag. For large M, this has substantial benefits, over the case that C\ is zero, but C 2 is the 
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cost associated with Algorithm [I] 

The specific method we use for generating a tag is given in Algorithm [2] We denote our graph as a set 
of vertices V and a set of edges E , and for any vefwe denote by E(v) the set of edges that have v as an 
endpoint. To each vertex v we assign a natural number i v . The central idea in Algorithm [2] is to search for 
the lexically lowest encoding of a representation of the connectivity matrix, where we are searching over a 
set of possible orderings of vertices. We used an MD5 hash of the connectivity matrix simply to use less 
memory. The requirement that the graph be planar is what allows us to generate the unique ordering of W 
at line twelve. We stored the configuration with the ordering of points that generated the lexically lowest 
encoding, i.e. the tag. 

When we generated a new stable configuration for a given N and s we would, when the new configuration 
had a non-degenerate extremal triangulation, also generate the associated tag. We would then collect the 
already generated configurations for that N and s whose energies were close to the energy of the newly 
generated configuration. Within the subset of these with unique extremal triangulations, we would search for 
the tag associated with the new configuration. If we found it, then because we stored the configurations with 
the orderings of points that generated the tag, we knew the rotation and reflection necessary that would be the 
isometry. It was our experience that, when there was an isometry, this method found it immediately. Further, 
it was our experience that almost all of the configurations had non-degenerate extremal triangulations. 

If we did not find matching tags, then we used Algorithm [I] to search for an isometry between the 
newly generated configuration and the existing configurations with similar energies. We only characterized 
a configuration as new for a given N and s if every configuration with that N and s had an energy that was 
sufficiently different to ensure that there was no isometry or that the application of Algorithm [I] did not find 
an isometry. The graph-isomorphism technique sped the process of finding isometries when they existed, 
but it was never used by itself to determine if a configuration was new or isometric to an existing one. 

There are two reasons why one should not rely exclusively on isomorphisms of non-degenerate extremal 
triangulations. The first is that it is possible, although we didn’t see this case, that the graph has a non-trivial 
automorphism, but that the associated mapping of points is not a self-isometry. Put another way, there may 
be two orderings of the points that lead to the same lexically minimal tag. An indication of this would 
be that, at line twenty-two of Algorithm [2j Tag was not None and T - Tag, but, that there is no isometry 
between the configurations that preserves the orderings of points. An extremely simple example is a triangle 
where no two sides have the same length. It has no self-isometries, but six graph-automorphisms. Also, 
there is the remote possibility for collisions in the MD5 algorithm. 

The more significant reason that graph isomorphisms alone are not sufficient to identify stable configu¬ 
rations is that we found two distinct configurations that both have non-degenerate extremal triangulations, 
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ALGORITHM 2. Generating a tag for a planar graph. 

1: Tag None. 

2: for v e V do 

3: for e e E(v) do 

4: for r € {Clockwise, Counterclockwise} do 

5: Reset all indices i v for v e V. 

6: i v <— 1. 

7: i w <— 2, where the edge e joins the points v and w. 

8: n <— 3. 

9: while there is a vertex that has not been indexed do 

10: v the vertex with the lowest index that has an unindexed neighbor. 

11: y <— jc’s neighbor with the lowest index. 

12: W <r- set of neighbors of v ordered by r and starting with y. 

13: for w eW do 

14: if w has not been indexed then 

15: i w n. 

16: n <— n + 1. 

17: end if 

18: end for 

19: end while 

20: P <— the connectivity matrix of the graph where the vertices are ordered by their indexing. 

21: T <— the MD5 cryptographic hash of P. 

22: if Tag = None or T < Tag then 

23: Tag <- T. 

24: end if 

25: end for 

26: end for 

27: end for 


but whose graphs were isomorphic. For N = 102 and s = 2, the configurations with the fourth and fifth 
lowest energy have non-degenerate extremal triangulation with isomorphic graphs. The dual graphs, i.e. the 
Voronoi cells, are shown in Figure[2] The difference in energy is substantially more than the estimated error 
in the energy sums. The energies are 5582.2331644897 and 5582.2332117851 respectively. Algorithm |T] 
did not identify an isometry. Further, out of thousands of computer trials, we reproduced the fourth lowest 
configuration 205 times and the fifth lowest configuration 100 times - these were not rare configurations. 
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FIG. 2. On the left is the configuration with the fourth lowest energy for N = 102 and s = 2, on the right is 
the configuration with the fifth lowest energy. They have the non-degenerate extremal triangulation with isomorphic 
graphs, but are distinct stable configurations. The dark (blue in the online version) cell has seven edges, two of which 
in the upper left and lower left are extremely short. The reader will notice that these two edges differ in the image on 
the left and the right. 


IV. UNSEEN SPECIES 


After months of running on many compute cores, we found that for most N between 120 and 180 the 
rate at which we discovered new stable configurations was still far from zero. We took this as an indication 
that more trials would result in more distinct stable configurations, and that we had not seen all of them. In 
response, we aimed to estimate the number of stable configurations that we didn’t see in our trials. Such an 
estimate cannot be made without additional assumptions, which we shall make clear as we proceed. There 
is some precedent for trying to estimate the number of unseen species. For example, Efron and Thisted 
estimated the number of words Shakespeare knew E ^l, although their approach is more sophisticated than 
ours. 

In broad terms our approach is as follows: We first compute the Good-Turing Frequency described 
below. This is an estimate for the combined probability of all the configurations we did not see. In addition 
this method produces estimated probabilities for the S stable configurations we did see, [p \, p 2 ,..., ps }• We 
assume that when these estimated probabilities are sorted in decreasing order the tail has a certain analytic 
form, i.e. that there is a p(ri) so that p n - p(ri) for large n. We obtain p(ri) from the data, and use it to 
compute how many more configurations we would need for the sum of the probabilities of those unseen 
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configurations to agree with the Good-Turing Frequency. 

The first assumption is that the number of stable configurations is finite. While this seems intuitively 
true, the function /(x) - x sin(l/x) for x ^ 0 and 0 for x = 0 has infinitely many local minima on the closed 
unit interval, indicating that a proof that there are finitely many stable configurations will depend on domain 
specific information. 

A second assumption is that, were we to use a different gradient based optimization technique, the es¬ 
timated probabilities for configurations we observed wouldn’t be so different as to change dramatically the 
estimates for the unseen species. While we have no proof, our instincts are that the basins of attraction 
for gradient descent methods all are qualitatively the same, and that the initial random configurations were 
sufficiently disordered so as not to be “nearer” a particular subset of stable configurations. Indeed, efforts 
to avoid the preponderance of stable configurations while searching for the global minimum has lead re¬ 
searchers away from purely gradient-based methods such as the work by Morris, Deavon and Hg MDH96 , and 
the work by Lakhab and BernoussI LEB13 . 

We now briefly summarize I. J. Good’s description of the Good-Turing estimateP 22 ^. Suppose we 
perform T trials where we observed some number of distinct species - stable configurations in our case. We 
let n r denote the number of species that we observed r times, and so 

oo 

J]n r r=T. (3) 

r= 1 

We then ask, for a species that we saw r times, what is a reasonable estimate for the fraction of the population 
that consists of that species? The most straightforward estimate, r/T , has the drawback that the sum of the 
fractions is one, i.e. this estimates assumes that there were no unseen species. This is almost certainly 
wrong in our case. If the likelihood of seeing each species is described by a binomial distribution, which 
is reasonable in our case, Good arrives at the following estimate 60053 (§2 Eq. 15) for the probability of a 
species that was observed r times 

r + 1 Sr+iOr+i) 

T + 1 S T (n r ) 

Here, in a slight abuse of notation, St (nr) indicates the expectation value for the number of species that we 
would expect to see r times in T trials. To be applicable, Good makes the following approximation 

T + 1 

£r+iOv+i) ~ n' r+ 1 , ~ 1, 

where n' r is the smoothed number of species seen r times. The need for smoothing arises because at large r, 
i.e. for species that occurred many times, the discreteness of the measurement n r within T trials becomes 
apparent. Gale and Sampson provide a method, which we used, for smoothing the data GS95 . Figures 1 and 
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2 in that publication make clear the need for, and effect of, smoothing. The result is the following estimate 
for the probability of a species (configuration) that occurred r times in T trials 60053 (Eqs. 2 and 2’) 


r+in'. 


Pr - 


r +1 


T n' 


The estimated probability of all species that occurred r times, denoted p r , is 


r + 1 


Pr - 


n r+l 


Summing these estimated probabilities over all observed speciesP 22 ^ (Eqs. 7,8) gives 

oo 1 oo 1 / oo \ 

Yj pr = rZ^l (r+1) " T Yj n 'r r ~ n '\ 

r= 1 r—\ \r= 1 / 

If the smoothing process is performed so that Eq. ([3]) holds with n r replaced with n', then the combined 
probability of all of our observed species is given by 

1 n\ 

-(T-n[) = 1--^ 

and so the estimate for the probability of the unseen species is 


P o = — • 


This method gives us estimated probabilities for each of the configurations we have observed. For 
N = 180 and s = 3 we show these estimated probabilities in Figure [3] The stair-step nature for the 
low estimated probabilities is an artifact of the finite number of samples. The estimated probability for a 
configuration depends only on the number of times the configuration occurred. This number can only be 
1,2,.... Gale and Sampson’s smoothing technique addresses a similar problem, in that the observed values 
of n r for large r must also be integral. We use a simple smoothing technique where, for a given probability, 
we take the geometric mean of the first and last configuration number as an estimate for the configuration 
number where that probability would occur. This is shown in Figure |4| 

We fit Ax + log b to the log-log tail of these data excluding the point corresponding to the configurations 
that occurred once. We are operating on two assumptions here: first, the tail of the probability distribution 
can be approximated p(n ) « bn A , and second, that the last point in Figure [ 4 ] is not “below” the fit line, as 
much as it is “to the left” of the fit line. That is to say, for this sample, we believe that many more unseen 
configurations whose probability is close to p\ than there are unseen configurations whose probability is 
close to or higher than P 2 . Another statement of this assumption is that nearly all the unseen configurations 
have probability less than p 2 , but not necessarily less than p\. 
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FIG. 3. This plot shows the estimated probability for all of the observed configurations, ordered in decreasing 
probability. 
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FIG. 4. This plot shows the smoothed estimated probability for all of the observed configurations. 
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FIG. 5. The number of observed distinct stable configurations and estimates for the total number of distinct stable 
configurations as a function of TV for s = 0. 


a 

o 

•a 

2 

a 

W) 

a 

a 

o 

o 


<D 

a 

£ 


100000 


10000 


1000 


100 


s=l observed 




jt 


•r 


■T 


** 

** ■ 

.****1 ** 

*_U 






II 

-f f l 

i 


*****' 


. 51 ^ 1010 ^ 0 * 0 ^ 


*********** 


10 


100 110 120 130 140 150 160 170 180 

Number of points N 


FIG. 6. The number of observed distinct stable configurations and estimates for the total number of distinct stable 
configurations as a function of TV for s = 1. 
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FIG. 7. The number of observed distinct stable configurations and estimates for the total number of distinct stable 
configurations as a function of TV for s = 2. 
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With this in mind we solve the following for Tf 

r T f 

Po + n\pi= I bn A dn , (4) 

JTi 

where po is the estimated combined probability of all the species we didn’t see, p\ is the estimated proba¬ 
bility for the species we saw once, n\ is the number of species that we saw once, and T; is the number of 
configurations we saw at least twice. Note that we are only guaranteed to get a value of Tf if A > -1. When 
we apply this method we get an estimate for the total number of stable configurations. These estimates as 
a function of N are plotted in Figures [5j [6|[7] and [8] for s = 0,1,2 and 3 respectively. In these figures we’ve 
only plotted results where the error was less then the value itself. 

If these estimates for the number of stable configurations are reasonable, and if the growth in the number 
of stable configurations is exponential, then fits from N - 100,..., 180 for the number of stable configura¬ 
tions as a function of N and s indicate that the number of stable configurations as a function of N and s is 
given by 


M(N,s) = C s N e ', 


where e s and C s are given by 


e 0 = 0.0741345 + 0.002804 
ei = 0.0789298 ± 0.00176 
<? 2 = 0.0836987 + 0.002186 
e 3 = 0.0878486 + 0.001698 


C 0 = exp(-3.91164 + 0.3044), 
Ci = exp(-3.97635 + 0.1992), 
C 2 = exp(-4.23711 ± 0.2583), 
C 3 = exp(-4.52585 + 0.1882). 


If we perform a similar fit to the number of observed configurations for s = 1, as opposed to the number 
of estimated configurations, we obtain M(N, 1) = (.31701 + .1) exp(.0518 + .0012 A), which is similar to 
Erber’s and Hockney’s estimate of M(N, 1) = 0.382 exp(0.0497A/^) noted above. This growth is considerably 
slower than the growth in the estimated number of stable configurations, which we feel is likely closer to 
the actual growth in the number of stable configurations. 


V. CONCLUSIONS 

When searching for isomorphisms between a set of configurations of points on a sphere, the use of a 
simple invariant under isometry, the discrete energy E s , quickly filtered out many configurations as not- 
isometric. After this a more comprehensive invariant under isometry, the graph of non-degenerate extremal 
triangulations, was extremely effective. 
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It is reasonable to express concern over the number of assumptions and over the sensitivity of Nf in 
Eq. ([4]) to the other parameters. The defensible conclusion is that there are substantially more stable config¬ 
urations than those we observed, and just as performing enough trials to observe the configuration with the 
lowest energy is a formidable technical challenge, so too is finding all the stable configurations. 
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